#ifndef VIRIAL_RADIUS_FINDER_H
#define VIRIAL_RADIUS_FINDER_H

#include <iostream>
#include <string>

#include "GadgetParticleType/gpt.01.h"
#include "Sorter/sorter.01.h"

class VRF {
	
	public:
		
		
		double M;
		double Epot;
		double Rvir;
		double Mvir;
		double Vvir;
		
		double G;
		double H0;
		double rho_crit;
		
		static const double pi = 3.141592654;		
		
		VRF() :
		G(43.0071),
		H0(100),
		rho_crit(3*H0*H0/(8*pi*G))
		{}
		
		// stellardynamicist's aproach 
		template <typename T1, typename T2>
		double find(uint N, T1 *m, T2 *pot){
			
			double Epot_2x(0);
			double M2(0); // squared total mass
			double m2(0); // sum of squared masses
			
			for (uint i=0; i<N; i++){
				Epot_2x+=pot[i];
				m2+=SQR(m[i]);
				M2+=m[i];
			}
			
			M = M2;
			M2 = SQR(M2);
			Epot = 0.5*Epot_2x;
			
			Rvir = -(M2-m2)/Epot_2x;
			

			return Rvir;
			
		}
		
		
		template <typename T>
		double find(uint N, P *p, T *pot){
			
			float *m = new float[N];
			
			for (uint i=0; i<N; i++)
			m[i] = p[i].Mass;
			
			find(N, m, pot);
			
			delete [] m;
		}
		
		
		// Cosmologist's aproach 
		template <typename T>
		double find(uint N, P *p, double delta, Vector<T> C){
			
			float *r = new float[N+100];
			uint *ind = new uint[N+100];
			
			
			for (uint i=0; i<N; i++) {
				r[i] = (*(FVector*)p[i].Pos-C).module();
				ind[i] = i;
			}
			
			Sorter sorter;
			sorter.sort(0, N-1, r, ind);
			
			double Mcum(0), Mdta(0);
			for (uint i=0; i<N; i++) {
				
				Mdta = delta*rho_crit*(4.0/3.0)*pi*r[i]*r[i]*r[i];
				Mcum+= p[ind[i]].Mass;
				
				if (Mdta-Mcum>0) { 
					Rvir = 0.5*(r[i]+r[i-1]); 
					Mvir = 0.5*(2*Mcum-p[ind[i]].Mass);
					Vvir = sqrt(G*Mvir/Rvir);  
					break; 
				}
				
			}
			
			delete [] r;
			delete [] ind;
			
			return Rvir;
			
		}
		
		void print(){
			
			printf("\nVirialRadiusFinder:\n\n");
			printf("Rvir = %g\n", Rvir);
			printf("Mvir = %g\n", Mvir);
			printf("Vvir = %g\n", Vvir);
			printf("rho_crit = %g\n", rho_crit);
			
		}


};


#endif